In this section we explored how the exposure and the infection by heat-stressed parasites affected the host reproduction. First we evaluated the number of individuals fully castrated in exposed uninfected and in infected Daphnias. This divisions into two sub analysis helped us to find more adapted models to describe the host reproductive patern.
One of the consequences of Pasteuria ramosa infection in Daphnia is the reduction in allocation to reproduction, which can also potentially be due to the energetic cost of immune response for uninfected hosts. Here we focused on the proportion of non-reproducing individuals.
FC1 = glmmTMB(data=infected, Never_repro ~ Isolate + Temperature + DeathAge, family = "binomial", REML=T)
FC1R = glmmTMB(data=infected, Never_repro ~ Isolate + Temperature + DeathAge + (1|Batch), family = "binomial", REML=T)
AICc(FC1, FC1R)
## df AICc
## FC1 7 239.8814
## FC1R 8 238.7983
FC1 = glmmTMB(data=infected, Never_repro ~ Isolate * Temperature + DeathAge, family = "binomial",
na.action = na.fail)
model_selection = dredge(FC1)
# see best models
model_selection
## Global model call: glmmTMB(formula = Never_repro ~ Isolate * Temperature + DeathAge,
## data = infected, family = "binomial", na.action = na.fail,
## ziformula = ~0, dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df logLik
## 16 0.79900 + -0.02228 + + + 13 -98.284
## 15 -0.19180 + + + + 12 -99.622
## 4 0.98420 + -0.01982 + 5 -108.183
## 3 0.11780 + + 4 -109.410
## 2 0.97250 + -0.02725 2 -111.825
## 8 0.94670 + -0.01938 + + 7 -107.700
## 7 0.09487 + + + 6 -108.859
## 1 -0.25280 + 1 -114.432
## 6 0.91880 + -0.02678 + 4 -111.381
## 5 -0.28880 + + 3 -113.875
## AICc delta weight
## 16 224.9 0.00 0.321
## 15 225.3 0.32 0.273
## 4 226.7 1.79 0.131
## 3 227.1 2.12 0.111
## 2 227.7 2.78 0.080
## 8 230.1 5.16 0.024
## 7 230.2 5.30 0.023
## 1 230.9 5.94 0.016
## 6 231.0 6.06 0.016
## 5 233.9 8.95 0.004
## Models ranked by AICc(x)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
| Never_repro | Never_repro | Never_repro | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p |
| DeathAge | 0.98 | 0.95 – 1.00 | 0.104 | 0.98 | 0.96 – 1.01 | 0.120 | |||
| Isolate [C1] | 2.27 | 0.66 – 7.83 | 0.193 | 2.46 | 0.72 – 8.38 | 0.151 | 1.23 | 0.49 – 3.10 | 0.658 |
| Isolate [C18] | 0.55 | 0.18 – 1.69 | 0.298 | 0.50 | 0.16 – 1.49 | 0.212 | 0.44 | 0.18 – 1.11 | 0.083 |
| Isolate [C24] | 0.65 | 0.20 – 2.09 | 0.472 | 0.66 | 0.21 – 2.08 | 0.474 | 0.50 | 0.20 – 1.28 | 0.150 |
| Temperature [linear] | 0.17 | 0.03 – 0.96 | 0.044 | 0.14 | 0.03 – 0.77 | 0.024 | |||
| Temperature [quadratic] | 0.28 | 0.07 – 1.12 | 0.073 | 0.29 | 0.07 – 1.16 | 0.080 | |||
|
Isolate [C1] × Temperature [linear] |
32.02 | 3.12 – 328.91 | 0.004 | 36.19 | 3.56 – 367.88 | 0.002 | |||
|
Isolate [C18] × Temperature [linear] |
2.04 | 0.26 – 16.14 | 0.497 | 2.80 | 0.37 – 21.15 | 0.319 | |||
|
Isolate [C24] × Temperature [linear] |
5.08 | 0.56 – 45.81 | 0.147 | 6.46 | 0.74 – 56.47 | 0.092 | |||
|
Isolate [C1] × Temperature [quadratic] |
10.98 | 1.60 – 75.39 | 0.015 | 10.14 | 1.50 – 68.66 | 0.018 | |||
|
Isolate [C18] × Temperature [quadratic] |
3.86 | 0.64 – 23.40 | 0.142 | 3.51 | 0.59 – 20.92 | 0.168 | |||
|
Isolate [C24] × Temperature [quadratic] |
3.68 | 0.59 – 23.04 | 0.163 | 3.31 | 0.54 – 20.28 | 0.196 | |||
| Observations | 167 | 167 | 167 | ||||||
| R2 Tjur | 0.178 | 0.163 | 0.074 | ||||||
best_model = get.models(model_selection, 2)[[1]]
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3977238 0.8536237 0.8579498 0.6620259 0.61213 0.9357635 0.05493602 0.2030844 0.1008657 0.9558345 0.1985479 0.8851658 0.4544661 0.1936561 0.1833452 0.2271874 0.5898748 0.8261575 0.5827213 0.7709696 ...
tab_model(best_model, show.intercept = F)
| Never_repro | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| Isolate [C1] | 2.46 | 0.72 – 8.38 | 0.151 |
| Isolate [C18] | 0.50 | 0.16 – 1.49 | 0.212 |
| Isolate [C24] | 0.66 | 0.21 – 2.08 | 0.474 |
| Temperature [linear] | 0.14 | 0.03 – 0.77 | 0.024 |
| Temperature [quadratic] | 0.29 | 0.07 – 1.16 | 0.080 |
|
Isolate [C1] × Temperature [linear] |
36.19 | 3.56 – 367.88 | 0.002 |
|
Isolate [C18] × Temperature [linear] |
2.80 | 0.37 – 21.15 | 0.319 |
|
Isolate [C24] × Temperature [linear] |
6.46 | 0.74 – 56.47 | 0.092 |
|
Isolate [C1] × Temperature [quadratic] |
10.14 | 1.50 – 68.66 | 0.018 |
|
Isolate [C18] × Temperature [quadratic] |
3.51 | 0.59 – 20.92 | 0.168 |
|
Isolate [C24] × Temperature [quadratic] |
3.31 | 0.54 – 20.28 | 0.196 |
| Observations | 167 | ||
| R2 Tjur | 0.163 | ||
## `summarise()` has grouped output by 'Temperature'. You can override using the
## `.groups` argument.
A heatwave experienced by the parasite appears to relax its castrating
effect on the host for most genotypes, with a strong opposit effect for
C1 genotyp, with a lower probability of total sterility in genotype C1,
which shows a dramatically increased risk of host sterility after
exposure to 40°C. This suggests an interaction between genotype and
heatwave experienced by the parasite.
options(na.action = "na.fail")
FC1 = glmmTMB(data=not_infected, Never_repro~ Isolate + Temperature + DeathAge, family = "binomial", REML=T)
FC1R = glmmTMB(data=not_infected, Never_repro ~ Isolate + Temperature+ DeathAge + (1|Batch), family = "binomial", REML=T)
AICc(FC1, FC1R)
## df AICc
## FC1 7 118.8705
## FC1R 8 106.9146
FC1 = glmmTMB(data=not_infected, Never_repro ~ Isolate * Temperature + DeathAge + (1|Batch), family = "binomial")
model_selection = dredge(FC1)
# Table of model comparison
model_selection
## Global model call: glmmTMB(formula = Never_repro ~ Isolate * Temperature + DeathAge +
## (1 | Batch), data = not_infected, family = "binomial", ziformula = ~0,
## dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df logLik
## 2 1.672 + -0.1281 3 -48.375
## 6 1.459 + -0.1208 + 5 -46.911
## 4 2.324 + -0.1529 + 6 -46.318
## 8 2.030 + -0.1347 + + 8 -45.341
## 16 13.600 + -0.8717 + + + 14 -40.057
## 5 -2.909 + + 4 -67.021
## 7 -2.961 + + + 7 -65.982
## 1 -2.830 + 2 -72.618
## 15 -3.165 + + + + 13 -61.318
## 3 -2.927 + + 5 -71.734
## AICc delta weight
## 2 102.9 0.00 0.499
## 6 104.1 1.24 0.268
## 4 105.0 2.17 0.168
## 8 107.4 4.51 0.052
## 16 110.2 7.32 0.013
## 5 142.2 39.37 0.000
## 7 146.5 43.63 0.000
## 1 149.3 46.43 0.000
## 15 150.4 47.56 0.000
## 3 153.8 50.89 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | Batch)
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3977238 0.8536237 0.8579498 0.6620259 0.61213 0.9357635 0.05493602 0.2030844 0.1008657 0.9558345 0.1985479 0.8851658 0.4544661 0.1936561 0.1833452 0.2271874 0.5898748 0.8261575 0.5827213 0.7709696 ...
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
| Never_repro | Never_repro | |||||
|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p |
| DeathAge | 0.88 | 0.82 – 0.94 | <0.001 | 0.89 | 0.83 – 0.95 | 0.001 |
| Temperature [linear] | 0.36 | 0.10 – 1.30 | 0.120 | |||
| Temperature [quadratic] | 0.61 | 0.17 – 2.17 | 0.442 | |||
| Random Effects | ||||||
| σ2 | 3.29 | 3.29 | ||||
| τ00 | 6.14 Batch | 5.52 Batch | ||||
| ICC | 0.65 | 0.63 | ||||
| N | 50 Batch | 50 Batch | ||||
| Observations | 218 | 218 | ||||
| Marginal R2 / Conditional R2 | 0.385 / 0.786 | 0.421 / 0.784 | ||||
# best model
best_model = get.models(model_selection, 1)[[1]]
summary(best_model)
## Family: binomial ( logit )
## Formula: Never_repro ~ DeathAge + (1 | Batch)
## Data: not_infected
##
## AIC BIC logLik deviance df.resid
## 102.8 112.9 -48.4 96.8 215
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Batch (Intercept) 6.137 2.477
## Number of obs: 218, groups: Batch, 50
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.67209 0.98865 1.691 0.09078 .
## DeathAge -0.12808 0.03582 -3.575 0.00035 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(best_model)
| Never_repro | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 5.32 | 0.77 – 36.96 | 0.091 |
| DeathAge | 0.88 | 0.82 – 0.94 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 Batch | 6.14 | ||
| ICC | 0.65 | ||
| N Batch | 50 | ||
| Observations | 218 | ||
| Marginal R2 / Conditional R2 | 0.385 / 0.786 | ||
Only the age at death seems to reprduce at least once in their life
time.
## `summarise()` has grouped output by 'DeathAgeGroup'. You can override using the
## `.groups` argument.
options(na.action = "na.fail")
dataC_ninf = subset(dataC, InfectedStatus!="infected")
summary(dataC_ninf)
## ID Temperature Isolate Clutch Nb_offspring
## Min. : 2 20: 84 C1 :45 Min. : 1.000 Min. : 0.00
## 1st Qu.: 218 30: 73 C14 :58 1st Qu.: 3.750 1st Qu.: 10.00
## Median : 668 40:100 C18 :67 Median : 8.000 Median : 47.00
## Mean : 588 C24 :60 Mean : 7.368 Mean : 43.02
## 3rd Qu.: 969 CTRL:27 3rd Qu.:11.000 3rd Qu.: 68.00
## Max. :1100 Max. :16.000 Max. :111.00
## NA's :29
## FirstClutchAge LastClutchAge DeathAge Accidental_death
## Min. : 1.00 Min. :10.00 Min. :15.00 Min. :1
## 1st Qu.:14.00 1st Qu.:39.00 1st Qu.:31.00 1st Qu.:1
## Median :18.00 Median :59.00 Median :62.00 Median :1
## Mean :21.77 Mean :50.47 Mean :51.47 Mean :1
## 3rd Qu.:27.00 3rd Qu.:67.00 3rd Qu.:67.00 3rd Qu.:1
## Max. :67.00 Max. :67.00 Max. :67.00 Max. :1
## NA's :30 NA's :32 NA's :250
## Volume_counted InfectedStatus Number of mature spores
## Min. :150 Length:257 Length:257
## 1st Qu.:150 Class :character Class :character
## Median :150 Mode :character Mode :character
## Mean :152
## 3rd Qu.:150
## Max. :200
## NA's :12
## Number of imature spores Nb of white spores cauliflowers grapes
## Length:257 Length:257 no :236 no :236
## Class :character Class :character NA's: 21 NA's: 21
## Mode :character Mode :character
##
##
##
##
## whites group statusDeath Number of immature spores
## Min. :0 C1440 : 29 Min. :0.0000 Length:257
## 1st Qu.:0 C2440 : 28 1st Qu.:0.0000 Class :character
## Median :0 NACTRL : 27 Median :1.0000 Mode :character
## Mean :0 C1840 : 26 Mean :0.5136
## 3rd Qu.:0 C1830 : 24 3rd Qu.:1.0000
## Max. :0 C130 : 19 Max. :1.0000
## (Other):104
## Never_repro Batch imm mat mature_load
## Min. :0.0000 2 : 8 Min. :0 Min. :0 Min. :0
## 1st Qu.:0.0000 17 : 8 1st Qu.:0 1st Qu.:0 1st Qu.:0
## Median :0.0000 28 : 8 Median :0 Median :0 Median :0
## Mean :0.1128 31 : 8 Mean :0 Mean :0 Mean :0
## 3rd Qu.:0.0000 3 : 7 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0
## Max. :1.0000 6 : 7 Max. :0 Max. :0 Max. :0
## (Other):211 NA's :12
## imm_load ratio repro_rate stade_max
## Min. :0 Min. : NA Min. :0.0000 None:257
## 1st Qu.:0 1st Qu.: NA 1st Qu.:0.3095
## Median :0 Median : NA Median :0.8358
## Mean :0 Mean :NaN Mean :0.7206
## 3rd Qu.:0 3rd Qu.: NA 3rd Qu.:1.0597
## Max. :0 Max. : NA Max. :1.6567
## NA's :12 NA's :257
FC1 = glmmTMB(data=dataC_ninf, Never_repro~ Isolate + Temperature + DeathAge, family = "binomial", REML=T)
FC1R = glmmTMB(data=dataC_ninf, Never_repro ~ Isolate + Temperature+ DeathAge + (1|Batch), family = "binomial", REML=T)
AICc(FC1, FC1R)
## df AICc
## FC1 8 134.3687
## FC1R 9 125.3196
FC1 = glmmTMB(data=not_infected, Never_repro ~ Isolate * Temperature + DeathAge + (1|Batch), family = "binomial")
model_selection = dredge(FC1)
# Table of model comparison
model_selection
## Global model call: glmmTMB(formula = Never_repro ~ Isolate * Temperature + DeathAge +
## (1 | Batch), data = not_infected, family = "binomial", ziformula = ~0,
## dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df logLik
## 2 1.672 + -0.1281 3 -48.375
## 6 1.459 + -0.1208 + 5 -46.911
## 4 2.324 + -0.1529 + 6 -46.318
## 8 2.030 + -0.1347 + + 8 -45.341
## 16 13.600 + -0.8717 + + + 14 -40.057
## 5 -2.909 + + 4 -67.021
## 7 -2.961 + + + 7 -65.982
## 1 -2.830 + 2 -72.618
## 15 -3.165 + + + + 13 -61.318
## 3 -2.927 + + 5 -71.734
## AICc delta weight
## 2 102.9 0.00 0.499
## 6 104.1 1.24 0.268
## 4 105.0 2.17 0.168
## 8 107.4 4.51 0.052
## 16 110.2 7.32 0.013
## 5 142.2 39.37 0.000
## 7 146.5 43.63 0.000
## 1 149.3 46.43 0.000
## 15 150.4 47.56 0.000
## 3 153.8 50.89 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | Batch)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
| Never_repro | Never_repro | |||||
|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p |
| DeathAge | 0.88 | 0.82 – 0.94 | <0.001 | 0.89 | 0.83 – 0.95 | 0.001 |
| Temperature [linear] | 0.36 | 0.10 – 1.30 | 0.120 | |||
| Temperature [quadratic] | 0.61 | 0.17 – 2.17 | 0.442 | |||
| Random Effects | ||||||
| σ2 | 3.29 | 3.29 | ||||
| τ00 | 6.14 Batch | 5.52 Batch | ||||
| ICC | 0.65 | 0.63 | ||||
| N | 50 Batch | 50 Batch | ||||
| Observations | 218 | 218 | ||||
| Marginal R2 / Conditional R2 | 0.385 / 0.786 | 0.421 / 0.784 | ||||
best_model = get.models(model_selection, 1)[[1]]
summary(best_model)
## Family: binomial ( logit )
## Formula: Never_repro ~ DeathAge + (1 | Batch)
## Data: not_infected
##
## AIC BIC logLik deviance df.resid
## 102.8 112.9 -48.4 96.8 215
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Batch (Intercept) 6.137 2.477
## Number of obs: 218, groups: Batch, 50
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.67209 0.98865 1.691 0.09078 .
## DeathAge -0.12808 0.03582 -3.575 0.00035 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(best_model)
| Never_repro | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 5.32 | 0.77 – 36.96 | 0.091 |
| DeathAge | 0.88 | 0.82 – 0.94 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 Batch | 6.14 | ||
| ICC | 0.65 | ||
| N Batch | 50 | ||
| Observations | 218 | ||
| Marginal R2 / Conditional R2 | 0.385 / 0.786 | ||
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.9088821 0.9008778 0.9317442 0.7483929 0.241452 0.2361889 0.7475164 0.4387352 0.390151 0.5425091 0.609537 0.7492443 0.6104253 0.7648734 0.5350931 0.2879339 0.134203 0.1612245 0.3157523 0.1570396 ...
Only the age at death seems to reprduce at least once in their life
time, the husbandry control inclusion does not change this.
## `summarise()` has grouped output by 'DeathAgeGroup'. You can override using the
## `.groups` argument.
In this section we were interested in understanding the role of parasite heat-stress and genotype on the number of offspring produced. Indeed, one of the consequences of infection in Daphnia is the reduction in allocation to parasite reproduction, (or a potential cost of immune response for uninfected individuals), even if the host is not fully castrated.
repro_inf = subset(data, data$Clutch>0 & InfectedStatus=="infected" & Temperature!="NA")
repro_inf$Temperature = factor(repro_inf$Temperature, ordered = T)
FC1 = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "nbinom2", REML=T)
#FC1 = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "poisson", REML=T)
#FC1 = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "gaussian", REML=T)
simulateResiduals(FC1, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.8741857 0.7996104 0.4441007 0.4213267 0.3856629 0.2460142 0.4371622 0.9950848 0.5756341 0.3335869 0.1840254 0.6387504 0.2376711 0.8736609 0.9228149 0.7387757 0.2361033 0.4922073 0.4105355 0.1956212 ...
FC1R = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge + (1|Batch), family = "nbinom2", REML=T)
AICc(FC1, FC1R)# no batch effect
## df AICc
## FC1 8 601.8980
## FC1R 9 601.5721
FC1 = glmmTMB(Nb_offspring ~ Isolate * Temperature + DeathAge ,
data = repro_inf, family = "nbinom2", na.action = na.fail)
model_selection = dredge(FC1)
model_selection
## Global model call: glmmTMB(formula = Nb_offspring ~ Isolate * Temperature + DeathAge,
## data = repro_inf, family = "nbinom2", na.action = na.fail,
## ziformula = ~0, dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df logLik
## 4 0.9519 + 0.03153 + 6 -283.958
## 2 0.2723 + 0.03680 3 -289.038
## 8 0.9455 + 0.03176 + + 8 -283.902
## 16 0.7038 + 0.02896 + + + 14 -277.332
## 6 0.3050 + 0.03628 + 5 -288.741
## 3 2.7120 + + 5 -291.710
## 7 2.7390 + + + 7 -291.682
## 15 2.2290 + + + + 13 -284.113
## 1 2.1350 + 2 -300.042
## 5 2.1400 + + 4 -299.296
## AICc delta weight
## 4 580.9 0.00 0.750
## 2 584.3 3.46 0.133
## 8 585.5 4.62 0.075
## 16 588.0 7.10 0.022
## 6 588.2 7.28 0.020
## 3 594.1 13.22 0.001
## 7 598.7 17.79 0.000
## 15 598.8 17.89 0.000
## 1 604.2 23.34 0.000
## 5 607.0 26.16 0.000
## Models ranked by AICc(x)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
| Nb_offspring | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| DeathAge | 1.03 | 1.02 – 1.05 | <0.001 |
| Isolate [C1] | 0.44 | 0.22 – 0.89 | 0.021 |
| Isolate [C18] | 0.48 | 0.28 – 0.82 | 0.008 |
| Isolate [C24] | 0.80 | 0.45 – 1.39 | 0.424 |
| Observations | 94 | ||
best_model = get.models(model_selection, 1)[[1]]
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.8445634 0.7757172 0.3188654 0.4678265 0.4046863 0.2574613 0.3273075 0.988557 0.6029993 0.2770807 0.2222652 0.5937299 0.200338 0.8887268 0.8925955 0.6897923 0.2627872 0.457771 0.3542096 0.1563993 ...
tab_model(best_model, show.intercept = F)
| Nb_offspring | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| DeathAge | 1.03 | 1.02 – 1.05 | <0.001 |
| Isolate [C1] | 0.44 | 0.22 – 0.89 | 0.021 |
| Isolate [C18] | 0.48 | 0.28 – 0.82 | 0.008 |
| Isolate [C24] | 0.80 | 0.45 – 1.39 | 0.424 |
| Observations | 94 | ||
summary(best_model)
## Family: nbinom2 ( log )
## Formula: Nb_offspring ~ DeathAge + Isolate + 1
## Data: repro_inf
##
## AIC BIC logLik deviance df.resid
## 579.9 595.2 -284.0 567.9 88
##
##
## Dispersion parameter for nbinom2 family (): 1.4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.951947 0.468983 2.030 0.04238 *
## DeathAge 0.031525 0.007622 4.136 3.54e-05 ***
## IsolateC1 -0.812114 0.352549 -2.304 0.02125 *
## IsolateC18 -0.738421 0.276490 -2.671 0.00757 **
## IsolateC24 -0.228773 0.286424 -0.799 0.42445
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
repro_nf = subset(data, data$Clutch>0 & InfectedStatus=="not_inf" & Temperature!="NA")
repro_nf$Temperature = factor(repro_nf$Temperature, ordered = T)
FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate * Temperature + DeathAge+ (1|Batch), family = "gaussian", REML=T)
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate * Temperature+ DeathAge + (1|Batch), family = "poisson")
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate * Temperature + DeathAge+ (1|Batch), family = "nbinom2", REML=T)
#FC1 = glmmTMB(data=repro_nf, log(Nb_offspring) ~ Isolate * Temperature + DeathAge + (1|Batch), family = "gaussian", REML=T)
FC1R = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate + Temperature + DeathAge + (1|Batch), family = "gaussian", REML=F)
FC1 = glmmTMB(Nb_offspring ~ Isolate + Temperature + DeathAge, data = repro_nf, family = "gaussian", REML=F)
AICc(FC1, FC1R)# Batch effect NOT important! But significantly improved the residuals
## df AICc
## FC1 8 1659.382
## FC1R 9 1660.230
FC1 = glmmTMB(Nb_offspring ~ Isolate * Temperature + DeathAge + (1|Batch), data = repro_nf, family = "gaussian", na.action = na.fail)
model_selection = dredge(FC1)
model_selection
## Global model call: glmmTMB(formula = Nb_offspring ~ Isolate * Temperature + DeathAge +
## (1 | Batch), data = repro_nf, family = "gaussian", na.action = na.fail,
## ziformula = ~0, dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df logLik
## 2 -22.63 + 1.306 4 -823.741
## 4 -20.37 + 1.313 + 7 -821.201
## 6 -22.45 + 1.294 + 6 -823.160
## 8 -20.41 + 1.303 + + 9 -820.623
## 5 48.57 + + 5 -913.919
## 1 50.26 + 3 -916.804
## 7 51.89 + + + 8 -913.084
## 3 53.77 + + 6 -915.627
## 15 52.47 + + + + 14 -912.123
## 16 -20.73 + 1.307 + + + 15
## AICc delta
## 2 1655.7 0.00
## 4 1657.0 1.31
## 6 1658.8 3.08
## 8 1660.2 4.54
## 5 1838.2 182.46
## 1 1839.7 184.04
## 7 1843.0 187.26
## 3 1843.7 188.01
## 15 1854.6 198.91
## 16
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | Batch)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
selected_models
## $`2`
## Formula: Nb_offspring ~ DeathAge + (1 | Batch)
## Data: repro_nf
## AIC BIC logLik df.resid
## 1655.4818 1668.5326 -823.7409 189
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## Batch (Intercept) 4.599
## Residual 16.721
##
## Number of obs: 193 / Conditional model: Batch, 50
##
## Dispersion estimate for gaussian family (sigma^2): 280
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) DeathAge
## -22.631 1.306
##
## $`4`
## Formula: Nb_offspring ~ DeathAge + Isolate + (1 | Batch)
## Data: repro_nf
## AIC BIC logLik df.resid
## 1656.4026 1679.2414 -821.2013 186
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## Batch (Intercept) 4.633
## Residual 16.482
##
## Number of obs: 193 / Conditional model: Batch, 50
##
## Dispersion estimate for gaussian family (sigma^2): 272
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) DeathAge IsolateC1 IsolateC18 IsolateC24
## -20.3723 1.3133 -4.2733 -0.6133 -6.5511
##
## $<NA>
## NULL
##
## attr(,"rank")
## function (x)
## do.call("rank", list(x))
## <environment: 0x137323810>
## attr(,"call")
## AICc(x)
## attr(,"class")
## [1] "function" "rankFunction"
## attr(,"beta")
## [1] "none"
best_model = get.models(model_selection, 1)[[1]]
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.78 0.656 0.584 0.568 0.468 0.564 0.596 0.448 0.512 0.388 0.464 0.368 0.504 0.452 0.672 0.536 0.708 0.592 0.356 0.488 ...
tab_model(best_model, show.intercept = F) # No effect
| Nb_offspring | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| DeathAge | 1.31 | 1.16 – 1.45 | <0.001 |
| Random Effects | |||
| σ2 | 279.59 | ||
| τ00 Batch | 21.15 | ||
| ICC | 0.07 | ||
| N Batch | 50 | ||
| Observations | 193 | ||
| Marginal R2 / Conditional R2 | 0.622 / 0.649 | ||
summary(best_model)
## Family: gaussian ( identity )
## Formula: Nb_offspring ~ DeathAge + (1 | Batch)
## Data: repro_nf
##
## AIC BIC logLik deviance df.resid
## 1655.5 1668.5 -823.7 1647.5 189
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Batch (Intercept) 21.15 4.599
## Residual 279.59 16.721
## Number of obs: 193, groups: Batch, 50
##
## Dispersion estimate for gaussian family (sigma^2): 280
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -22.63107 4.34365 -5.21 1.89e-07 ***
## DeathAge 1.30586 0.07379 17.70 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculation of daily reproduction rate
repro_nf$daily_repro_rate = repro_nf$Nb_offspring / repro_nf$DeathAge
mean(repro_nf$daily_repro_rate, na.rm = TRUE)
## [1] 0.8219201
(sem = sd(repro_nf$daily_repro_rate, na.rm = TRUE) / sqrt(length(repro_nf$daily_repro_rate)))
## [1] 0.02774078
Only the effect of time is important here; each day, a Daphnia produces on average 0.30 more offspring.
dataC$Isolate = as.character(dataC$Isolate)
dataC$Isolate[is.na(dataC$Isolate)] = "unexposed"
dataC$Isolate = as.factor(dataC$Isolate)
dataC$Temperature = as.character(dataC$Temperature)
dataC$Temperature[dataC$Temperature == "CTRL"] = 20
dataC$Temperature = factor(dataC$Temperature, ordered = TRUE)
repro_nfC = subset(dataC, dataC$Clutch>0 & InfectedStatus!="inf")
repro_nfC$Temperature = factor(repro_nfC$Temperature, labels = c("20", "30","40"), ordered = T)
FC1 = glmmTMB(data=repro_nfC, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "gaussian", REML=F)
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate + Temperature+ DeathAge, family = "poisson", REML=F)
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "nbinom2", REML=T)
simulateResiduals(FC1, plot=T) # I don't find better and I have no reason to exclude the outlier.
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.724 0.608 0.58 0.488 0.6 0.48 0.656 0.54 0.672 0.492 0.484 0.408 0.4 0.548 0.464 0.352 0.8 0.624 0.824 0.668 ...
FC1R = glmmTMB(data=repro_nfC, Nb_offspring ~ Isolate + Temperature + DeathAge + (1|Batch), family = "gaussian", REML=F)
FC1 = glmmTMB(Nb_offspring ~ Isolate + Temperature + DeathAge, data = repro_nfC, family = "gaussian", REML=F)
AICc(FC1, FC1R)# Batch effect NOT important and do not improve residuals anyway.
## df AICc
## FC1 9 1945.748
## FC1R 10 1946.820
FC1 = glmmTMB(Nb_offspring ~ Isolate * Temperature + DeathAge, data = repro_nfC, family = "gaussian", na.action = na.fail)
model_selection = dredge(FC1)
model_selection
## Global model call: glmmTMB(formula = Nb_offspring ~ Isolate * Temperature + DeathAge,
## data = repro_nfC, family = "gaussian", na.action = na.fail,
## ziformula = ~0, dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df logLik
## 2 -22.74 + 1.305 3 -966.679
## 4 -24.62 + 1.311 + 7 -963.982
## 6 -22.39 + 1.295 + 5 -966.198
## 8 -24.57 + 1.305 + + 9 -963.461
## 16 -24.28 + 1.314 + + + 15 -962.257
## 5 47.93 + + 4 -1078.562
## 1 48.49 + 2 -1083.107
## 7 50.74 + + + 8 -1077.345
## 3 52.70 + + 6 -1079.491
## 15 46.80 + + + + 14 -1076.437
## AICc delta weight
## 2 1939.5 0.00 0.682
## 4 1942.5 3.01 0.151
## 6 1942.7 3.20 0.137
## 8 1945.7 6.28 0.029
## 16 1956.8 17.31 0.000
## 5 2165.3 225.84 0.000
## 1 2170.3 230.80 0.000
## 7 2171.3 231.88 0.000
## 3 2171.4 231.90 0.000
## 15 2182.8 243.38 0.000
## Models ranked by AICc(x)
best_model = get.models(model_selection, 1)[[1]]
simulateResiduals(best_model, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.776 0.628 0.604 0.564 0.528 0.548 0.648 0.632 0.548 0.436 0.528 0.432 0.496 0.436 0.432 0.452 0.724 0.516 0.784 0.664 ...
selected_models = get.models(model_selection, 1:nrow(model_selection))
tab_model(selected_models, show.intercept = F) # No effect
| Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | Nb_offspring | |||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p |
| DeathAge | 1.30 | 1.18 – 1.43 | <0.001 | 1.31 | 1.18 – 1.44 | <0.001 | 1.29 | 1.17 – 1.42 | <0.001 | 1.31 | 1.18 – 1.43 | <0.001 | 1.31 | 1.18 – 1.44 | <0.001 | |||||||||||||||
| Isolate [C14] | 4.03 | -2.97 – 11.04 | 0.259 | 3.83 | -3.22 – 10.88 | 0.287 | 3.21 | -4.95 – 11.37 | 0.441 | 1.19 | -10.42 – 12.80 | 0.841 | 0.82 | -10.79 – 12.43 | 0.890 | 5.39 | -8.07 – 18.85 | 0.433 | ||||||||||||
| Isolate [C18] | 3.72 | -3.09 – 10.52 | 0.284 | 3.82 | -3.00 – 10.63 | 0.272 | 3.17 | -4.64 – 10.98 | 0.426 | -3.36 | -14.52 – 7.80 | 0.555 | -4.35 | -15.56 – 6.87 | 0.447 | 0.98 | -11.90 – 13.86 | 0.882 | ||||||||||||
| Isolate [C24] | -2.42 | -9.45 – 4.62 | 0.501 | -2.51 | -9.64 – 4.63 | 0.491 | -2.96 | -11.08 – 5.16 | 0.475 | -5.37 | -17.12 – 6.37 | 0.370 | -6.47 | -18.12 – 5.19 | 0.277 | -1.23 | -14.62 – 12.17 | 0.858 | ||||||||||||
| Isolate [CTRL] | 1.21 | -7.28 – 9.69 | 0.780 | 2.25 | -7.51 – 12.00 | 0.652 | -0.81 | -18.17 – 16.54 | 0.927 | -7.95 | -23.93 – 8.03 | 0.330 | -15.55 | -29.36 – -1.74 | 0.027 | 7.65 | -20.95 – 36.26 | 0.600 | ||||||||||||
| Temperature [linear] | 1.63 | -2.10 – 5.37 | 0.391 | 1.71 | -2.54 – 5.96 | 0.431 | -1.60 | -14.45 – 11.26 | 0.808 | 9.22 | 3.23 – 15.20 | 0.003 | 7.35 | 0.41 – 14.29 | 0.038 | 19.74 | -1.19 – 40.67 | 0.065 | ||||||||||||
| Temperature [quadratic] | 0.81 | -3.19 – 4.80 | 0.692 | 0.98 | -3.23 – 5.18 | 0.649 | 1.01 | -8.88 – 10.90 | 0.841 | -2.19 | -8.71 – 4.33 | 0.510 | -1.09 | -8.01 – 5.83 | 0.757 | -8.20 | -24.46 – 8.05 | 0.323 | ||||||||||||
|
Isolate [C14] × Temperature [linear] |
3.67 | -11.49 – 18.83 | 0.635 | -14.13 | -38.99 – 10.72 | 0.265 | ||||||||||||||||||||||||
|
Isolate [C18] × Temperature [linear] |
3.79 | -11.12 – 18.70 | 0.618 | -14.78 | -39.20 – 9.63 | 0.235 | ||||||||||||||||||||||||
|
Isolate [C24] × Temperature [linear] |
3.84 | -10.96 – 18.64 | 0.611 | -12.61 | -36.88 – 11.67 | 0.309 | ||||||||||||||||||||||||
|
Isolate [C14] × Temperature [quadratic] |
-2.63 | -15.81 – 10.55 | 0.696 | 8.53 | -13.14 – 30.20 | 0.440 | ||||||||||||||||||||||||
|
Isolate [C18] × Temperature [quadratic] |
4.00 | -8.12 – 16.11 | 0.518 | 9.27 | -10.71 – 29.24 | 0.363 | ||||||||||||||||||||||||
|
Isolate [C24] × Temperature [quadratic] |
-2.94 | -16.34 – 10.46 | 0.667 | 5.90 | -16.16 – 27.96 | 0.600 | ||||||||||||||||||||||||
| Observations | 228 | 228 | 228 | 228 | 228 | 228 | 228 | 228 | 228 | 228 | ||||||||||||||||||||
| R2 / R2 adjusted | 0.640 / 0.637 | 0.648 / 0.639 | 0.641 / 0.635 | 0.650 / 0.637 | 0.654 / 0.631 | 0.039 / 0.026 | 0.000 / -0.004 | 0.049 / 0.019 | 0.031 / 0.009 | 0.057 / -0.000 | ||||||||||||||||||||
summary(best_model)
## Family: gaussian ( identity )
## Formula: Nb_offspring ~ DeathAge + 1
## Data: repro_nfC
##
## AIC BIC logLik deviance df.resid
## 1939.4 1949.6 -966.7 1933.4 225
##
##
## Dispersion estimate for gaussian family (sigma^2): 282
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -22.74119 3.70948 -6.131 8.76e-10 ***
## DeathAge 1.30473 0.06482 20.127 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Again, the husbandry control do not change the conclusions.
Here we just described the average number of offspring produced per day in total in each group
# Function to compute mean and 95% confidence interval
mean_ci = function(x) {
x = x[!is.na(x)]
m = mean(x)
se = sd(x) / sqrt(length(x))
c(mean = m,
CI_low = m - 1.96 * se,
CI_high = m + 1.96 * se)
}
# Define groups
g1 = mean_ci(repro_nfC$repro_rate[repro_nfC$Isolate == "CTRL"])
g2 = mean_ci(repro_nf$repro_rate[repro_nfC$Temperature == 20])
g3 = mean_ci(repro_nf$repro_rate)
g4 = mean_ci(repro_inf$repro_rate)
rbind(
noninfected_unexposed = g1,
noninfected_exposed_20 = g2,
noninfected_exposed_all = g3,
infected_all = g4
)
## mean CI_low CI_high
## noninfected_unexposed 0.7203259 0.5747725 0.8658793
## noninfected_exposed_20 0.7463513 0.6490932 0.8436094
## noninfected_exposed_all 0.8219201 0.7675482 0.8762921
## infected_all 0.1699311 0.1302886 0.2095736